Code
# Load required packages
library(plm)
library(fixest)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(patchwork)In my previous post on Fixed vs Between vs Random Effects, I explained three different approaches to panel data analysis. The core strength of fixed effects models is their ability to “control for all time-invariant characteristics of counties.” However, this approach addresses only one dimension of potential confounding in panel data—it misses the crucial time dimension.
Two-way fixed effects (TWFE) models extend one-way fixed effects by controlling for both entity-specific and time-specific unobserved factors. This dual control makes TWFE particularly powerful for causal inference in policy evaluation.
While one-way fixed effects control for time-invariant characteristics of counties, two-way fixed effects simultaneously control for:
\[Y_{it} = \beta X_{it} + \alpha_i + \lambda_t + \epsilon_{it}\]
Where:
County Fixed Effects Control For:
Time Fixed Effects Control For:
# Load required packages
library(plm)
library(fixest)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(patchwork)Let’s create a simple county-level dataset with clear county and time effects, using 12 counties observed over 8 years (2015-2022).
Dependent variable: Percentage of uninsured population in each county
Independent variable: Medicaid expansion indicator (binary treatment variable)
Control variable: County unemployment rate
County fixed effects: Represent persistent differences between counties due to factors like:
Time fixed effects: Capture year-specific shocks affecting all counties equally, such as:
This setup allows us to demonstrate how TWFE isolates the causal effect of Medicaid expansion by controlling for both stable county characteristics and common time trends that could otherwise confound our estimates.
set.seed(123)
n_counties <- 12
n_years <- 8
total_obs <- n_counties * n_years
# Generate identifiers
county_id <- rep(1:n_counties, each = n_years)
year <- rep(2015:2022, times = n_counties)
# Simple county fixed effects (different baseline levels)
county_effects <- rep(rnorm(n_counties, mean = 0, sd = 20), each = n_years)
# Simple time fixed effects (common trends affecting all counties)
time_effects_values <- c(-5, -3, 0, 2, 3, 1, -2, 4) # One for each year
time_effects <- rep(time_effects_values, times = n_counties)
# Treatment: Medicaid expansion (some counties expand in 2018, others never)
treatment_year <- rep(sample(c(2018, NA), n_counties,
replace = TRUE, prob = c(0.5, 0.5)), each = n_years)
medicaid_expansion <- ifelse(is.na(treatment_year), 0,
ifelse(year >= treatment_year, 1, 0))
# Control variable
unemployment <- rep(rnorm(n_counties, mean = 5, sd = 1), each = n_years) +
rnorm(total_obs, 0, 0.5)
# Outcome: Uninsured rate
treatment_effect <- -4 # True effect: 4 percentage point reduction
uninsured_rate <- 12 + # Baseline
treatment_effect * medicaid_expansion + # Treatment effect
0.3 * unemployment + # Unemployment effect
county_effects + # County differences
time_effects + # Time trends
rnorm(total_obs, 0, 1) # Noise
# Create dataset
twfe_data <- data.frame(
county_id = factor(county_id),
year = year,
medicaid_expansion = medicaid_expansion,
uninsured_rate = uninsured_rate,
unemployment = unemployment,
treatment_year = treatment_year
)
# Show sample
head(twfe_data, 12) %>%
kable(digits = 2, caption = "Sample Two-Way Fixed Effects Data") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| county_id | year | medicaid_expansion | uninsured_rate | unemployment | treatment_year |
|---|---|---|---|---|---|
| 1 | 2015 | 0 | -2.20 | 5.91 | 2018 |
| 1 | 2016 | 0 | -0.47 | 5.55 | 2018 |
| 1 | 2017 | 0 | 1.67 | 6.15 | 2018 |
| 1 | 2018 | 1 | 0.56 | 6.14 | 2018 |
| 1 | 2019 | 1 | 3.07 | 6.11 | 2018 |
| 1 | 2020 | 1 | 0.06 | 6.05 | 2018 |
| 1 | 2021 | 1 | -3.37 | 5.98 | 2018 |
| 1 | 2022 | 1 | 2.07 | 5.67 | 2018 |
| 2 | 2015 | 0 | 1.66 | 4.37 | 2018 |
| 2 | 2016 | 0 | 6.83 | 4.34 | 2018 |
| 2 | 2017 | 0 | 7.19 | 4.18 | 2018 |
| 2 | 2018 | 1 | 7.46 | 4.42 | 2018 |
The visualization shows three key patterns:
# Add treatment group labels
twfe_data <- twfe_data %>%
mutate(treatment_group = ifelse(is.na(treatment_year), "Never Treated", "Treated 2018"))
# Plot 1: Raw trends show both county differences and time patterns
ggplot(twfe_data, aes(x = year, y = uninsured_rate, color = treatment_group)) +
geom_line(aes(group = county_id), alpha = 0.5) +
stat_summary(fun = mean, geom = "line", size = 2) +
geom_vline(xintercept = 2018, linetype = "dashed", color = "black") +
scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
labs(title = "Raw Data: County Trends Over Time",
subtitle = "Shows county differences, time trends, and treatment effects mixed together",
x = "Year", y = "Uninsured Rate (%)", color = "Group") +
theme_minimal()# Plot 2: Average by year (shows time effects)
time_means <- twfe_data %>%
group_by(year) %>%
summarize(avg_uninsured = mean(uninsured_rate))
ggplot(time_means, aes(x = year, y = avg_uninsured)) +
geom_line(size = 1.5, color = "darkgreen") +
geom_point(size = 3, color = "darkgreen") +
labs(title = "Average Uninsured Rate by Year",
subtitle = "Shows common time trends affecting all counties",
x = "Year", y = "Average Uninsured Rate (%)") +
theme_minimal()# Plot 3: Average by county (shows county effects)
county_means <- twfe_data %>%
group_by(county_id) %>%
summarize(avg_uninsured = mean(uninsured_rate)) %>%
mutate(county_rank = rank(avg_uninsured))
ggplot(county_means, aes(x = county_rank, y = avg_uninsured)) +
geom_col(fill = "steelblue", alpha = 0.7) +
labs(title = "Average Uninsured Rate by County",
subtitle = "Shows persistent differences between counties",
x = "County (Ranked by Uninsured Rate)", y = "Average Uninsured Rate (%)") +
theme_minimal()The key to understanding TWFE is seeing how it removes both county and time effects through “demeaning”:
# Step 1: Calculate means
twfe_data <- twfe_data %>%
group_by(county_id) %>%
mutate(county_mean = mean(uninsured_rate)) %>%
ungroup() %>%
group_by(year) %>%
mutate(year_mean = mean(uninsured_rate)) %>%
ungroup() %>%
mutate(
overall_mean = mean(uninsured_rate),
# Step 2: Apply two-way demeaning transformation
uninsured_demeaned = uninsured_rate - county_mean - year_mean + overall_mean
)The two-way fixed effects transformation follows this formula:
\(Y_{it}^* = Y_{it} - \bar{Y}_i - \bar{Y}_t + \bar{Y}\)
Where:
Step 1: Why subtract the county mean (\(\bar{Y}_i\))?
# Visualization of each step
ggplot(twfe_data, aes(x = year, y = uninsured_rate, color = treatment_group)) +
geom_line(aes(group = county_id), alpha = 0.6) +
stat_summary(fun = mean, geom = "line", size = 2) +
scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
labs(title = "Step 1: Original Data",
subtitle = "Raw uninsured rates with all variation",
x = "Year", y = "Uninsured Rate (%)", color = "Group") +
theme_minimal()Each county has its own baseline level due to permanent characteristics (geography, demographics, institutions). By subtracting each county’s average, we remove these time-invariant differences. After this step, all counties are centered around zero, eliminating the county fixed effects.
Step 2: Why subtract the year mean (\(\bar{Y}_t\))?
# After removing county effects
twfe_data <- twfe_data %>%
mutate(after_county_demean = uninsured_rate - county_mean + overall_mean)
ggplot(twfe_data, aes(x = year, y = after_county_demean, color = treatment_group)) +
geom_line(aes(group = county_id), alpha = 0.6) +
stat_summary(fun = mean, geom = "line", size = 2) +
scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
labs(title = "Step 2: After Removing County Effects",
subtitle = "County differences eliminated, time patterns remain",
x = "Year", y = "County-Demeaned Rate", color = "Group") +
theme_minimal()Now we can see that the distance between county lines has narrowed considerably! After removing county-specific effects, the data reveals a clear common time pattern across all counties.
This pattern reflects year-specific factors that affect all counties equally - such as economic cycles, federal healthcare policies, national insurance market trends, and macro-level health crises. By subtracting each year’s average from the county-demeaned data, we remove these time-specific shocks, eliminating the systematic year-to-year variation that would otherwise confound our treatment effect estimates.
What remains after this second demeaning step is variation that differs across counties within the same time period - this is the key identifying variation that allows TWFE to isolate causal effects while controlling for both persistent county differences and common temporal trends.
Step 3: Why add back the overall mean (\(\bar{Y}\))?
# After removing both county and time effects
ggplot(twfe_data, aes(x = year, y = uninsured_demeaned, color = treatment_group)) +
geom_line(aes(group = county_id), alpha = 0.6) +
stat_summary(fun = mean, geom = "line", size = 2) +
scale_color_manual(values = c("Never Treated" = "blue3", "Treated 2018" = "red3")) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.7) +
labs(title = "Step 3: After Removing County & Time Effects",
subtitle = "Only treatment effects and idiosyncratic variation remain",
x = "Year", y = "Two-Way Demeaned Rate", color = "Group") +
theme_minimal()Now we see a clearer pattern emerge between never-treated counties and treated counties! The treatment effect becomes much more apparent after removing both county-specific baselines and common time trends.
This transformation works through careful mathematical adjustment: when we subtract both county means (\(\bar{Y}_i\)) and year means (\(\bar{Y}_t\)), we inadvertently subtract the overall mean (\(\bar{Y}\)) twice. This double-subtraction occurs because the overall mean is embedded within both county averages and year averages. Adding back \(\bar{Y}\) corrects this mathematical artifact, ensuring our transformation removes only the county-specific and time-specific components while preserving the interpretable scale of the original data.
The result is data that isolates the treatment effect by comparing how treated counties differ from control counties at the same point in time, after accounting for both persistent county characteristics and common temporal shocks.
After this transformation, only two sources of variation remain:
This is the “identifying variation” - variation that occurs within counties and within time periods simultaneously. TWFE uses this remaining variation to estimate treatment effects while controlling for all county-specific and time-specific confounders.
Consider county A with high baseline uninsured rates and county B with low baseline rates. In a recession year, both counties see increases in uninsured rates.
This is why TWFE provides such powerful control for confounding factors.
Again, two-way fixed effects estimation follows a sequential demeaning process that differs from conventional one-way fixed effects:
Step 1: Individual demeaning (Same as conventional FE) \[\tilde{Y}_{it} = Y_{it} - \bar{Y}_i\] \[\tilde{X}_{it} = X_{it} - \bar{X}_i\]
This removes county-specific means from both the outcome and explanatory variables.
In our example: We subtract each county’s average uninsured rate and average Medicaid expansion status from their respective yearly values. This eliminates persistent differences between counties - some counties that are always high-uninsured vs. low-uninsured due to demographics, healthcare infrastructure, or political culture.
Step 2: Time demeaning of demeaned variables (Major difference from conventional FE) \[\ddot{Y}_{it} = \tilde{Y}_{it} - \bar{\tilde{Y}}_t\] \[\ddot{X}_{it} = \tilde{X}_{it} - \bar{\tilde{X}}_t\]
This removes time-specific means from the already county-demeaned variables. Note that \(\bar{\tilde{Y}}_t\) is the average of the county-demeaned outcomes in year \(t\).
In our example: We subtract the yearly averages of the county-demeaned data. This removes common time trends affecting all counties - such as economic cycles, federal policy changes, or national healthcare trends that cause all counties’ uninsured rates to move together over time.
Step 3: OLS on double-demeaned variables \[\ddot{Y}_{it} = \beta \ddot{X}_{it} + \ddot{\epsilon}_{it}\]
Run standard OLS regression on the double-demeaned variables. The coefficient \(\beta\) gives us the two-way fixed effects estimate.
In our example: We estimate how changes in Medicaid expansion (relative to county and time averages) affect changes in uninsured rates (relative to county and time averages). This isolates the pure treatment effect by comparing treated counties to control counties at the same point in time, after removing both county-specific baselines and time-specific shocks.
The key insight is that we don’t just subtract raw year means - we subtract the means of the already county-demeaned data. This ensures we’re removing time effects from data that has already had county effects removed.
In our context: This means we’re not just controlling for the fact that 2020 was a bad year for healthcare (raw time effect), but rather how 2020 affected counties after accounting for their baseline differences. This prevents time effects from being confounded with county-level characteristics.
Mathematically, this is equivalent to the single-step formula we saw earlier: \[Y_{it}^* = Y_{it} - \bar{Y}_i - \bar{Y}_t + \bar{Y}\]
But the two-step process makes it clearer how TWFE builds on conventional fixed effects methodology.
Now let’s estimate all four specifications using the fixest package, which provides fast and modern fixed effects estimation:
# Model 1: No fixed effects (pooled OLS)
model_none <- feols(uninsured_rate ~ medicaid_expansion + unemployment,
data = twfe_data)
# Model 2: County fixed effects only
model_county <- feols(uninsured_rate ~ medicaid_expansion + unemployment | county_id,
data = twfe_data)
# Model 3: Time fixed effects only
model_time <- feols(uninsured_rate ~ medicaid_expansion + unemployment | year,
data = twfe_data)
# Model 4: Two-way fixed effects (county + time)
model_twfe <- feols(uninsured_rate ~ medicaid_expansion + unemployment | county_id + year,
data = twfe_data)
# Display all models in a comparison table
etable(model_none, model_county, model_time, model_twfe,
headers = c("No FE", "County FE", "Time FE", "Two-Way FE"),
digits = 3,
se.below = TRUE) model_none model_county model_time model_twfe
No FE County FE Time FE Two-Way FE
Dependent Var.: uninsured_rate uninsured_rate uninsured_rate uninsured_rate
Constant 44.4***
(8.69)
medicaid_expansion -11.5** 0.712* -23.0*** -3.52***
(3.53) (0.273) (0.727) (0.334)
unemployment -5.13** 0.104 -5.19*** -0.176
(1.82) (0.736) (0.846) (0.166)
Fixed-Effects: -------------- -------------- -------------- --------------
county_id No Yes No Yes
year No No Yes Yes
__________________ ______________ ______________ ______________ ______________
S.E. type IID by: county_id by: year by: county_id
Observations 96 96 96 96
R2 0.17989 0.98135 0.31740 0.99753
Within R2 -- 0.01276 0.30731 0.42858
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The fixest package uses a clean syntax where fixed effects are specified after the | symbol:
y ~ x | fe1: One-way fixed effects for fe1y ~ x | fe1 + fe2: Two-way fixed effects for fe1 and fe2y ~ x: No fixed effects (standard OLS)# Extract and compare treatment effects specifically
models_list <- list(
"No FE" = model_none,
"County FE" = model_county,
"Time FE" = model_time,
"Two-Way FE" = model_twfe
)
# Create coefficient comparison dataframe
coef_comparison <- data.frame(
Model = names(models_list),
Coefficient = sapply(models_list, function(x) coef(x)["medicaid_expansion"]),
SE = sapply(models_list, function(x) se(x)["medicaid_expansion"])
) %>%
mutate(
CI_Lower = Coefficient - 1.96 * SE,
CI_Upper = Coefficient + 1.96 * SE,
True_Effect = -4, # Our known true effect
Bias = round(Coefficient - True_Effect, 3)
)
# Display detailed comparison
coef_comparison %>%
mutate(across(c(Coefficient, SE, CI_Lower, CI_Upper), ~round(.x, 3))) %>%
kable(caption = "Treatment Effect Estimates: Medicaid Expansion on Uninsured Rate") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Model | Coefficient | SE | CI_Lower | CI_Upper | True_Effect | Bias | |
|---|---|---|---|---|---|---|---|
| No FE.medicaid_expansion | No FE | -11.452 | 3.529 | -18.368 | -4.535 | -4 | -7.452 |
| County FE.medicaid_expansion | County FE | 0.712 | 0.273 | 0.176 | 1.248 | -4 | 4.712 |
| Time FE.medicaid_expansion | Time FE | -22.995 | 0.727 | -24.419 | -21.570 | -4 | -18.995 |
| Two-Way FE.medicaid_expansion | Two-Way FE | -3.523 | 0.334 | -4.177 | -2.869 | -4 | 0.477 |
# Visualize the comparison with specified order
coef_comparison$Model <- factor(coef_comparison$Model,
levels = c("No FE", "County FE", "Time FE", "Two-Way FE"))
ggplot(coef_comparison, aes(x = Model, y = Coefficient)) +
geom_hline(yintercept = -4, color = "grey40", linetype = "dashed", size = 1.2) +
geom_point(size = 4, color = "red", shape = 15) +
geom_errorbar(aes(ymin = CI_Lower, ymax = CI_Upper),
width = 0.2, size = 1, color = "black") +
coord_flip() +
scale_x_discrete(limits = rev(levels(coef_comparison$Model))) + # Reverse for coord_flip
labs(title = "Treatment Effect Estimates Across Specifications",
subtitle = "Grey dashed line shows true simulated effect (-4 percentage points)",
x = "Model Specification",
y = "Estimated Effect on Uninsured Rate (pp)",
caption = "Error bars show 95% confidence intervals") +
theme_minimal() +
theme(plot.title = element_text(face = "bold", size = 14))This visualization shows how different model specifications dramatically affect our estimates:
No Fixed Effects: Severely overestimates the effect (-12 pp) due to uncontrolled selection bias - counties choosing Medicaid expansion differ systematically from non-expanding counties.
County FE Only: Underestimates (-2 pp) by ignoring time trends. Fails to account for broader healthcare improvements happening during the study period.
Time FE Only: Dramatically overestimates (-22 pp) by ignoring county differences. Counties with better healthcare systems were more likely to expand, confounding results.
Two-Way FE: Provides the most credible estimate (-4 pp), closest to our true simulated effect, by controlling for both county characteristics and time trends simultaneously.
The fixest package allows us to easily extract and examine the estimated fixed effects:
county_fe <- fixef(model_twfe)$county_id
# Display county fixed effects
county_fe_df <- data.frame(
County_ID = names(county_fe),
County_FE = round(county_fe, 2)
) %>%
arrange(County_FE)
ggplot(county_fe_df, aes(x = reorder(County_ID, County_FE), y = County_FE)) +
geom_col(fill = "steelblue", alpha = 0.7) +
coord_flip() +
labs(title = "County Fixed Effects",
subtitle = "Persistent differences between counties (relative to reference)",
x = "County ID", y = "Fixed Effect Estimate") +
theme_minimal()County Fixed Effects: These represent persistent differences between counties after controlling for time trends and other covariates.
time_fe <- fixef(model_twfe)$year
# Display time fixed effects
time_fe_df <- data.frame(
Year = names(time_fe),
Time_FE = round(time_fe, 2)
)
ggplot(time_fe_df, aes(x = Year, y = Time_FE)) +
geom_col(fill = "darkgreen", alpha = 0.7) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.7) +
labs(title = "Time Fixed Effects",
subtitle = "Year-specific shocks affecting all counties",
x = "Year", y = "Fixed Effect Estimate") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))Time Fixed Effects: These capture year-specific shocks that affected all counties equally.
Key Insight: By estimating and removing these fixed effects, our treatment effect estimate is based only on variation that differs across counties within the same time period - providing cleaner causal identification.
The results show how different specifications can lead to different conclusions:
The two-way fixed effects estimate should be closest to our true simulated effect of -4 percentage points.